#include <stdio.h>
#include <ctype.h>
#include <stdlib.h>
#include <syslog.h>
#include <string.h>
#include <float.h>
#include <math.h>
#include "phylip.h"
#include "dist.h"
#include "fitch.h"

#define MAXNUMTREES   100000000 /* a number bigger than conceivable numtrees */
#define epsilonf        0.000001   /* a very small but not too small number  */
#define delta           0.0001      /* a not quite so small number */
#define zsmoothings     10    /* number of zero-branch correction iterations */

FILE *infile;
FILE *outtree;
char infilename[FILE_NAME_LEN];
char outtreename[FILE_NAME_LEN];
/** minimum evolution */
boolean minev;
/** global rearrangements */
boolean global;
/** randomize input order of species */
boolean jumble;
/** number of times to jumble */
boolean njumble;
/** use lengths from user trees */
boolean lengths;
/** lower-triangular data matrix */
boolean lower;
/** upper-triangular data matrix */
boolean upper;
/** negative branch lengths allowed */
boolean negallowed;
/** outgroup root */
boolean outgropt;
/** subreplicates */
boolean replicates;
/** outgroup species number */
int outgrno;
int inseed;
int nums;
int jumb = 0;
longer seed;
double power;
int *enterorder;
int nonodes2;
intvector *reps;
vector *x;
node *addwhere;
tree curtree, priortree, bestree, bestree2;
/**
 * From each of those places go back to all others 
 * nx comes from firsttraverse 
 * sum comes from evaluate via firsttraverse 
 */
void secondtraverse(node *q, double y, int *nx, double *sum)
{
	double z=0.0, TEMP=0.0;

	z = y + q->v;
  	if (q->tip) 
	{
    	TEMP = q->d[(*nx) - 1] - z;
    	*sum += q->w[(*nx) - 1] * (TEMP * TEMP);
  	} 
	else 
	{
    	secondtraverse(q->next->back, z, nx, sum);
    	secondtraverse(q->next->next->back, z, nx,sum);
  	}
}
/**
 * Go through tree calculating branch lengths 
 */
void firsttraverse(node *p, int *nx, double *sum)
{
  	if (minev && (p != curtree.start))
    	*sum += p->v;
  	if (p->tip) 
	{
    	if (!minev) 
		{
      		*nx = p->index;
      		secondtraverse(p->back, 0.0, nx, sum);
      	}
  	} 
	else 
	{
    	firsttraverse(p->next->back, nx,sum);
    	firsttraverse(p->next->next->back, nx,sum);
  	}
} 
/**
 * Evaluate likelihood of a tree 
 */
double evaluate(tree *t)
{
  	double sum=0.0;
  	int nx=0;
  
  	firsttraverse(t->start->back, &nx, &sum);
  	firsttraverse(t->start, &nx, &sum);
  	if ((!minev) && replicates && (lower || upper))
    	sum /= 2;
  	t->likelihood = -sum;
  	return (-sum);
}
/**
 * Compute distance between an interior node and tips 
 */
void nudists(node *x, node *y)
{
	int nq=0, nr=0, nx=0, ny=0;
  	double dil=0, djl=0, wil=0, wjl=0, vi=0, vj=0;
  	node *qprime, *rprime;

  	qprime = x->next;
  	rprime = qprime->next->back;
  	qprime = qprime->back;
  	ny = y->index;
  	dil = qprime->d[ny - 1];
  	djl = rprime->d[ny - 1];
  	wil = qprime->w[ny - 1];
  	wjl = rprime->w[ny - 1];
  	vi = qprime->v;
  	vj = rprime->v;
  	x->w[ny - 1] = wil + wjl;
  	if (wil + wjl <= 0.0)
    	x->d[ny - 1] = 0.0;
  	else
    	x->d[ny - 1] = ((dil - vi) * wil + (djl - vj) * wjl) / (wil + wjl);
  	nx = x->index;
  	nq = qprime->index;
  	nr = rprime->index;
  	dil = y->d[nq - 1];
  	djl = y->d[nr - 1];
  	wil = y->w[nq - 1];
  	wjl = y->w[nr - 1];
  	y->w[nx - 1] = wil + wjl;
  	if (wil + wjl <= 0.0)
    	y->d[nx - 1] = 0.0;
  	else
    	y->d[nx - 1] = ((dil - vi) * wil + (djl - vj) * wjl) / (wil + wjl);
}
void filltraverse(node *pb, node *qb, boolean contin)
{
  	if (qb->tip)
    	return;
  	if (contin) 
	{
		filltraverse(pb, qb->next->back,contin);
		filltraverse(pb, qb->next->next->back,contin);
		nudists(qb, pb);
		return;
  	}
  	if (!qb->next->back->tip)
    	nudists(qb->next->back, pb);
  	if (!qb->next->next->back->tip)
    	nudists(qb->next->next->back, pb);
}
void fillin(node *pa, node *qa, boolean contin)
{
  	if (!pa->tip) 
	{
    	fillin(pa->next->back, qa, contin);
    	fillin(pa->next->next->back, qa, contin);
  	}
  	filltraverse(pa, qa, contin);
}
/**
 * Traverse updating these views 
 */
void alter(node *x, node *y)
{
  	nudists(x, y);
  	if (!y->tip) 
	{
    	alter(x, y->next->back);
    	alter(x, y->next->next->back);
  	}
}
/**
 * Renew information about subtrees 
 */
void nuview(node *p)
{
  	int i=0;
  	node *q, *r, *pprime, *temp;

  	q = p->next;
  	r = q->next;
  	for (i = 1; i <= 3; i++) 
	{
    	temp = p;
    	pprime = p->back;
    	alter(p, pprime);
    	p = q;
    	q = r;
    	r = temp;
  	}
}
/** 
 * Compute distances among three neighbors of a node 
 */
void makedists(node *p)
{
  	int i=0, nr=0, ns=0;
  	node *q, *r, *s;

  	r = p->back;
  	nr = r->index;
  	for (i = 1; i <= 3; i++) 
	{
		q = p->next;
		s = q->back;
		ns = s->index;
		if (s->w[nr - 1] + r->w[ns - 1] <= 0.0)
			p->dist = 0.0;
		else
			p->dist = (s->w[nr - 1] * s->d[nr - 1] + r->w[ns - 1] * r->d[ns - 1]) /
				(s->w[nr - 1] + r->w[ns - 1]);
		p = q;
		r = s;
		nr = ns;
  	}
}
/** 
 * Make new branch length 
 */
void makebigv(node *p)
{
  	int i=0;
  	node *temp, *q, *r;

  	q = p->next;
  	r = q->next;
  	for (i = 1; i <= 3; i++) 
	{
    	if (p->iter) 
		{
      		p->v = (p->dist + r->dist - q->dist) / 2.0;
      		p->back->v = p->v;
    	}
    	temp = p;
    	p = q;
    	q = r;
    	r = temp;
  	}
} 
/** 
 * Iterate branch lengths if some are to be zero 
 */
void correctv(node *p)
{
  	node *q, *r, *temp;
  	int i=0, j=0, n=0, nq=0, nr=0, ntemp=0;
  	double wq=0.0, wr=0.0;

  	q = p->next;
  	r = q->next;
  	n = p->back->index;
  	nq = q->back->index;
  	nr = r->back->index;
  	for (i = 1; i <= zsmoothings; i++) 
	{
    	for (j = 1; j <= 3; j++) 
		{
      		if (p->iter) 
			{
        		wr = r->back->w[n - 1] + p->back->w[nr - 1];
        		wq = q->back->w[n - 1] + p->back->w[nq - 1];
        		if (wr + wq <= 0.0 && !negallowed)
          			p->v = 0.0;
        		else
          			p->v = ((p->dist - q->v) * wq + (r->dist - r->v) * wr) / (wr + wq);
        		if (p->v < 0 && !negallowed)
          			p->v = 0.0;
        		p->back->v = p->v;
      		}
			temp = p;
			p = q;
			q = r;
			r = temp;
			ntemp = n;
			n = nq;
			nq = nr;
			nr = ntemp;
    	}
  	}
}
/** 
 * Update branch lengths around a node 
 */
void update(node *p)
{
	if (p->tip)
    	return;
  	makedists(p);
  	if (p->iter || p->next->iter || p->next->next->iter) 
	{
    	makebigv(p);
    	correctv(p);
  	}
  	nuview(p);
}
/** 
 * Go through tree getting new branch lengths and views 
 */
void smooth(node *p)
{
  	if (p->tip)
    	return;
  	update(p);
  	smooth(p->next->back);
  	smooth(p->next->next->back);
}
/**
 * Put p and q together and iterate info. on resulting tree 
 */
void insert_(node *p, node *q, boolean contin_)
{
	double x=0.0, oldlike;
  	hookup(p->next->next, q->back);
  	hookup(p->next, q);
  	x = q->v / 2.0;
  	p->v = 0.0;
  	p->back->v = 0.0;
  	p->next->v = x;
  	p->next->back->v = x;
  	p->next->next->back->v = x;
  	p->next->next->v = x;
  	fillin(p->back, p, contin_);
  	evaluate(&curtree);
  	do 
	{
    	oldlike = curtree.likelihood;
    	smooth(p);
    	smooth(p->back);
    	evaluate(&curtree);
  	} 
	while (fabs(curtree.likelihood - oldlike) > delta);
}
/**
 * Initialize branch lengths and views in a tip 
 */
void setuptipf(int m, tree *t)
{
	int i=0;
  	intvector n=(int *)malloc(spp * sizeof(int)); 
  	node *WITH;
	if ( n == NULL )
	{
		syslog(LOG_ERR,"setuptipf: malloc failed\n");
		exit( -1 );
	}
  	WITH = t->nodep[m - 1];
  	memcpy(WITH->d, x[m - 1], (nonodes2 * sizeof(double)));
 	memcpy(n, reps[m - 1], (spp * sizeof(int)));
  	for (i = 0; i < spp; i++) 
	{
    	if (i + 1 != m && n[i] > 0) 
		{
      		if (WITH->d[i] < epsilonf)
        	WITH->d[i] = epsilonf;
      		WITH->w[i] = n[i] / exp(power * log(WITH->d[i]));
    	}
		else 
		{
      		WITH->w[i] = 0.0;
     		WITH->d[i] = 0.0;
    	}
  	}
  	for (i = spp; i < nonodes2; i++) 
	{
    	WITH->w[i] = 1.0;
    	WITH->d[i] = 0.0;
  	}
  	WITH->index = m;
  	if (WITH->iter) WITH->v = 0.0;
  		free(n);
}
/**
 * Initialize and hook up a new tip 
 */
void buildnewtip(int m, tree *t, int nextsp)
{
  	node *p;
  	setuptipf(m, t);
  	p = t->nodep[nextsp + spp - 3];
  	hookup(t->nodep[m - 1], p);
} 
/**
 * Make and initialize a three-species tree 
 */
void buildsimpletree(tree *t, long nextsp)
{
  	curtree.start=curtree.nodep[enterorder[0] - 1]; 
  	setuptipf(enterorder[0], t);
  	setuptipf(enterorder[1], t);
  	hookup(t->nodep[enterorder[0] - 1], t->nodep[enterorder[1] - 1]);
  	buildnewtip(enterorder[2], t, nextsp);
  	insert_(t->nodep[enterorder[2] - 1]->back, t->nodep[enterorder[0] - 1],
		false);
}
/** read in species name (was in phylip.c) */
void initname(int i)
{
  	int j;
  	for (j = 0; j < nmlngth; j++) 
	{
		if (eoff(infile) | eoln(infile))
		{
		  	syslog(LOG_ERR,"end-of-line or end-of-file"
		  	" in the middle of species name for species %d\n", i+1);
		  	exit(-1);
		}
		nayme[i][j] = gettc(infile);
		if ((nayme[i][j] == '(') || (nayme[i][j] == ')') || (nayme[i][j] == ':')
		    || (nayme[i][j] == ',') || (nayme[i][j] == ';') || (nayme[i][j] == '[')
		    || (nayme[i][j] == ']')) {
		  	syslog(LOG_ERR,"Species name may not contain characters ( ) : ; , [ ] \n"
		  		"       In name of species number %d there is character %c\n",
		    	i+1, nayme[i][j]);
		  	exit(-1);
		}
	}
}
/** Fitch: read species number (was in phylip.c) */
void inputnumbers2(int *spp, int *nonodes, int n)
{
  	if (fscanf(infile, "%d", spp) != 1 || *spp <= 0) 
	{
   		syslog(LOG_ERR,"Unable to read the number of species in data set\n"
    	"The input file is incorrect (perhaps it was not saved text only).\n");
  	}
  	*nonodes = *spp * 2 - n;
}
/**
 * Allocate remaining arrays etc.
 */
void allocrest()
{
  	int i;
	x = (vector *)malloc(spp*sizeof(vector));
  	if ( x == NULL )
		goto bail;
	reps = (intvector *)malloc(spp*sizeof(intvector));
  	if ( reps == NULL )
		goto bail;
	for (i=0;i<spp;++i)
	{
		x[i]=(vector)malloc(nonodes2 * sizeof(double));
		if ( x[i] == NULL )
			goto bail;
		reps[i]=(intvector)malloc(spp * sizeof(int));
		if ( reps[i] == NULL )
			goto bail;
  	}
  	nayme = (naym *)malloc(spp*sizeof(naym));
	if ( nayme == NULL )
		goto bail;
  	enterorder = (int *)malloc(spp*sizeof(int));
	if ( enterorder == NULL )
		goto bail;
	for (i=0;i<spp;++i)
		enterorder[i]=0;
	return;
	bail:
		syslog(LOG_ERR,"allocrest: failed to allocate memory\n");
		exit( -1 );
}
/** initialize variables */
static void doinit()
{
	inputnumbers2(&spp, &nonodes2, 1);
	alloctree(&curtree.nodep, nonodes2);
	allocd(nonodes2, curtree.nodep);
	allocw(nonodes2, curtree.nodep);
	alloctree(&bestree.nodep, nonodes2);
	allocd(nonodes2, bestree.nodep);
	allocw(nonodes2, bestree.nodep);
	alloctree(&priortree.nodep, nonodes2);
	allocd(nonodes2, priortree.nodep);
	allocw(nonodes2, priortree.nodep);
	if (njumble > 1) 
	{
		alloctree(&bestree2.nodep, nonodes2);
		allocd(nonodes2, bestree2.nodep);
		allocw(nonodes2, bestree2.nodep);
	}
  	allocrest();
}
/**
 * Parse the argument list and set globals
 * Currently: [-i infile] [-o outtree]
 * The following booleans are defined. Values are Y or N. 
 * For definitions see code comments.
 * [-d minev] [-g global] [-j jumble] [-n lengths]
 * [-l lower] [-s replicates] [-o outgropt]
 * [-u upper] [-a negallowed] 
 * Floating point args: [-p power]
 * No arguments are required. If unspecified, the 
 * infile is stdin and the outtree is stdout. 
 * @return 1 if the args were OK, 0 otherwise
 */
static int readArgs( int argc, char ** argv )
{
	int i;
	int inseed0=0;
	for ( i=1;i<argc;i++ )
	{
		if ( strlen(argv[i]) > 1 )
		{
			if ( argv[i][0] == '-' && argc > i+1 )
			{
				switch ( argv[i][1] )
				{
					case 'i':
						strncpy(infilename,argv[i+1],FILE_NAME_LEN);
						break;
					case 'd':
						minev = (toupper(argv[i+1][0])=='Y')?true:false;
						break;
					case 'g':
						global = (toupper(argv[i+1][0])=='Y')?true:false;
						break;
					case 'j':
						jumble = (toupper(argv[i+1][0])=='Y')?true:false;
						if (jumble)
            				initseed(&inseed, &inseed0, seed);
					break;
					case 'n':
						lengths = (toupper(argv[i+1][0])=='Y')?true:false;
						break;
					case 'l':
						lower = (toupper(argv[i+1][0])=='Y')?true:false;
						break;
					case 'r':
						upper = (toupper(argv[i+1][0])=='Y')?true:false;
						break;
					case 's':
						replicates = (toupper(argv[i+1][0])=='Y')?true:false;
						break;
					case 'a':
						negallowed = (toupper(argv[i+1][0])=='Y')?true:false;
						break;
					case 'o':
						outgropt = (toupper(argv[i+1][0])=='Y')?true:false;
						break;
					case 'p':
						sscanf(argv[i+1],"%lf",&power);
						break;
					default:
						return 0;
						break;
				}
			}
		}
	}
	return (upper&&lower)?0:1;
}
/**
 * Set sensible defaults for all the global options
 */
static void setOptionDefaults()
{
	minev = false;
	global = false;
	jumble = false;
	njumble = 1;
	lengths = false;
	lower = false;
	upper = false;
	negallowed = false;
	outgropt = false;
	replicates = false;
	outgrno = 1;
	power = 2.0;
}
/** 
 * Make a copy of a node 
 */
void copynode(node *c, node *d)
{
  	memcpy(d->d, c->d, nonodes2*sizeof(double));
  	memcpy(d->w, c->w, nonodes2*sizeof(double));
  	d->v = c->v;
  	d->iter = c->iter;
  	d->dist = c->dist;
  	d->xcoord = c->xcoord;
  	d->ycoord = c->ycoord;
  	d->ymin = c->ymin;
  	d->ymax = c->ymax;
}
/** 
 * Make copy of a tree a to tree b 
 */
static void copy_(tree *a, tree *b)
{
	int i, j=0;
  	node *p, *q;

  	for (i = 0; i < spp; i++) 
	{
    	copynode(a->nodep[i], b->nodep[i]);
    	if (a->nodep[i]->back) 
		{
      		if (a->nodep[i]->back == a->nodep[a->nodep[i]->back->index - 1])
        		b->nodep[i]->back = b->nodep[a->nodep[i]->back->index - 1];
      		else if (a->nodep[i]->back
            	== a->nodep[a->nodep[i]->back->index - 1]->next)
        	b->nodep[i]->back = b->nodep[a->nodep[i]->back->index - 1]->next;
      		else
        		b->nodep[i]->back
          		= b->nodep[a->nodep[i]->back->index - 1]->next->next;
    	}
    	else b->nodep[i]->back = NULL;
  	}
  	for (i = spp; i < nonodes2; i++) 
	{
		p = a->nodep[i];
		q = b->nodep[i];
    	for (j = 1; j <= 3; j++) 
		{
      		copynode(p, q);
      		if (p->back) 
			{
        		if (p->back == a->nodep[p->back->index - 1])
          			q->back = b->nodep[p->back->index - 1];
        		else if (p->back == a->nodep[p->back->index - 1]->next)
          			q->back = b->nodep[p->back->index - 1]->next;
        		else
          			q->back = b->nodep[p->back->index - 1]->next->next;
      		}
      		else
        		q->back = NULL;
      		p = p->next;
      		q = q->next;
    	}
  	}
  	b->likelihood = a->likelihood;
  	b->start = a->start;
}
/** 
 * Traverse through a tree, finding best place to add p 
 */
void addtraverse(node *p, node *q, boolean contin, int *numtrees,
	boolean *succeeded)
{
	insert_(p, q, true);
  	(*numtrees)++;
  	if (evaluate(&curtree) > (bestree.likelihood + 
		epsilonf * fabs(bestree.likelihood)))
	{
    	copy_(&curtree, &bestree);
    	addwhere = q;
    	(*succeeded)=true;
  	}
  	copy_(&priortree, &curtree);
  	if (!q->tip && contin) 
	{
    	addtraverse(p, q->next->back, contin,numtrees,succeeded);
    	addtraverse(p, q->next->next->back, contin,numtrees,succeeded);
  	}
}
/**
 * Re_move p and record in q where it was 
 */
void re_move(node **p, node **q)
{
  	*q = (*p)->next->back;
  	hookup(*q, (*p)->next->next->back);
  	(*p)->next->back = NULL;
  	(*p)->next->next->back = NULL;
  	update(*q);
  	update((*q)->back);
}
/**
 * Does global rearrangements 
 */
void globrearrange(int* numtrees,boolean* succeeded) 
{
  	tree globtree;
  	tree oldtree;
  	int i,j,k,num_sibs,num_sibs2;
  	node *where,*sib_ptr,*sib_ptr2;
  	double oldbestyet = curtree.likelihood;
  	int success = false;
 
  	alloctree(&globtree.nodep,nonodes2);
  	alloctree(&oldtree.nodep,nonodes2);
  	setuptree(&globtree,nonodes2);
  	setuptree(&oldtree,nonodes2);
  	allocd(nonodes2, globtree.nodep);
  	allocd(nonodes2, oldtree.nodep);
  	allocw(nonodes2, globtree.nodep);
  	allocw(nonodes2, oldtree.nodep);
  	copy_(&curtree,&globtree);
  	copy_(&curtree,&oldtree);
  	for ( i = spp ; i < nonodes2 ; i++ ) 
	{
    	num_sibs = count_sibs(curtree.nodep[i]);
    	sib_ptr  = curtree.nodep[i];
		/* this appears to be mere progress-reporting, 
		   especially as stdout is not treeout */
    	/*if ( (i - spp) % (( nonodes2 / 72 ) + 1 ) == 0 )
      		putchar('.');
    	fflush(stdout);*/
    	for ( j = 0 ; j <= num_sibs ; j++ ) 
		{
      		re_move(&sib_ptr,&where);
      		copy_(&curtree,&priortree);
      
      		if (where->tip) 
			{
        		copy_(&oldtree,&curtree);
        		copy_(&oldtree,&bestree);
        		sib_ptr=sib_ptr->next;
        		continue;
      		}
      		else num_sibs2 = count_sibs(where);
      		sib_ptr2 = where;
      		for ( k = 0 ; k < num_sibs2 ; k++ ) 
			{
        		addwhere = NULL;
        		addtraverse(sib_ptr,sib_ptr2->back,true,numtrees,succeeded);
        		if ( addwhere && where != addwhere && where->back != addwhere
              		&& bestree.likelihood > globtree.likelihood) 
				{
            		copy_(&bestree,&globtree);
            		success = true;
        		}
        		sib_ptr2 = sib_ptr2->next;
      		} 
      		copy_(&oldtree,&curtree);
      		copy_(&oldtree,&bestree);
      		sib_ptr = sib_ptr->next;
    	}
  	}
  	copy_(&globtree,&curtree);
  	copy_(&globtree,&bestree);
  	if (success && globtree.likelihood > oldbestyet)  
	{
    	*succeeded = true;
  	}
  	else  
	{
    	*succeeded = false;
  	}
  	freed(nonodes2, globtree.nodep);
  	freed(nonodes2, oldtree.nodep);
  	freew(nonodes2, globtree.nodep);
  	freew(nonodes2, oldtree.nodep);
  	freetree(&globtree.nodep,nonodes2);
  	freetree(&oldtree.nodep,nonodes2);
}
void rearrange(node *p, int *numtrees, int *nextsp, boolean *succeeded)
{
  	node *q, *r;
  	if (!p->tip && !p->back->tip) 
	{
    	r = p->next->next;
    	re_move(&r, &q);
    	copy_(&curtree, &priortree);
    	addtraverse(r, q->next->back, false, numtrees,succeeded);
    	addtraverse(r, q->next->next->back, false, numtrees,succeeded);
    	copy_(&bestree, &curtree);
		/*
    	if (global && ((*nextsp) == spp)) 
		{
      		putchar('.');
      		fflush(stdout);
    	}
		*/
  	}
  	if (!p->tip) 
	{
    	rearrange(p->next->back, numtrees,nextsp,succeeded);
    	rearrange(p->next->next->back, numtrees,nextsp,succeeded);
  	}
} 
/** construct the tree */
void maketree()
{
	int nextsp,numtrees;
  	boolean succeeded=false;
  	int i;
	if (jumb == 1) 
	{
		inputdata(replicates, lower, upper, x, reps);
		setuptree(&curtree, nonodes2);
		setuptree(&priortree, nonodes2);
		setuptree(&bestree, nonodes2);
		if (njumble > 1) 
			setuptree(&bestree2, nonodes2);
	}
	for (i = 1; i <= spp; i++)
		enterorder[i - 1] = i;
	if (jumble)
		randumize(seed, enterorder);
	nextsp = 3;
	buildsimpletree(&curtree, nextsp);
	curtree.start = curtree.nodep[enterorder[0] - 1]->back;
	if (jumb == 1) 
		numtrees = 1;
	nextsp = 4;
	while (nextsp <= spp) 
	{
		nums = nextsp;
		buildnewtip(enterorder[nextsp - 1], &curtree, nextsp);
		copy_(&curtree, &priortree);
		bestree.likelihood = -DBL_MAX;
		curtree.start = curtree.nodep[enterorder[0] - 1]->back;
		addtraverse(curtree.nodep[enterorder[nextsp - 1] - 1]->back,
			curtree.start, true, &numtrees,&succeeded);
		copy_(&bestree, &curtree);
		succeeded = true;
		while (succeeded) 
		{
			succeeded = false;
			curtree.start = curtree.nodep[enterorder[0] - 1]->back;
			if (nextsp == spp  && global)
  				globrearrange (&numtrees,&succeeded);
			else
			{
  				rearrange(curtree.start,&numtrees,&nextsp,&succeeded);
			}
		}
		if (njumble > 1) 
		{
			if (jumb == 1 && nextsp == spp)
  				copy_(&bestree, &bestree2);
			else if (nextsp == spp) 
			{
  				if (bestree2.likelihood < bestree.likelihood)
    				copy_(&bestree, &bestree2);
			}
		}
      	if (nextsp == spp && jumb == njumble)
		{
        	int col = 0;
        	if (njumble > 1)
				copy_(&bestree2, &curtree);
        	curtree.start = curtree.nodep[outgrno - 1]->back;
    		treeout(curtree.start, &col, 0.43429445222, true,
				curtree.start);
      	}
		nextsp++;
	}
} 
/**
 * Main entry point. Completely rewritten to remove dependencies 
 * on X11 etc and to remove interactive commandline interface. 
 * Now we read from stdin and write to stdout by default, making 
 * it easy to pipe stuff through the program. We also eliminate the 
 * outfile and only write an outtree.
 */
int main(int argc, char *argv[])
{
	setOptionDefaults();
	if ( readArgs(argc,argv) )
	{
		if ( strlen(infilename) > 0 )
			infile = fopen( infilename, "r" );
		else 
			infile = stdin;
		if ( infile != NULL )
		{
			if ( strlen(outtreename) > 0 )
				outtree = fopen( outtreename, "w" );
			else
				outtree = stdout;
			if ( outtree != NULL )
			{
  				doinit();
				for (jumb = 1; jumb <= njumble; jumb++)
    				maketree();
				if (eoln(infile))
  					scan_eoln(infile);
				fclose( outtree );
			}
			else
				syslog(LOG_ERR,"Couldn't open file %s for writing\n",
					outtreename );
			fclose( infile );
		}
		else
			syslog(LOG_ERR,"Couldn't open file %s for reading\n",
				infilename );
	}
	else
		syslog(LOG_ERR,"Drawtree: invalid argument list\n");
	
  	return 0;
}
